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The  (2-5-2)  Spline  Function 


Jae  H.  Park  and  Leonard  A.  Ferrari 


Abstract.  Splines  have  been  used  extensively  in  the  interpolation  of  mul- 
tidimensional data  sets.  Linear  interpolation  utilizes  second  order  splines 
(first  degree  piecewise  polynomials)  and  has  widespread  popularity  be- 
cause of  its  ease  of  implementation.  Cubic  splines  are  often  used  when 
higher  degrees  of  smoothness  are  required  of  the  interpolation  process. 
Linear  interpolation  has  the  advantages  of  not  requiring  the  solution  of 
an  inverse  problem  (the  data  points  are  themselves  the  coefficients  of  the 
triangular  basis  functions)  and  extremely  efficient  generation  of  the  output 
sample  points.  Unfortunately,  the  linear-interpolating  function  has  only 
C°  continuity  (the  function  is  continuous  but  its  derivatives  are  discontin- 
uous) and  therefore  lacks  the  required  smoothness  for  many  applications. 
We  provide  a new  algorithm  in  this  paper  based  on  the  efficient  derivative 
summation  approach  to  spline  rendering.  Cubic  B-spline  interpolation  for 
uniformly  spaced  data  points  provides  C2  continuity.  The  interpolation 
function  can  be  rendered  quite  efficiently  from  the  basis  coefficients  and 
the  basis  function,  using  a cascade  of  four  running  average  filters.  Unser 
et  al.  have  shown  a digital  filter  solution  for  the  inverse  problem  of  ob- 
taining the  spline  coefficients  from  the  data  points.  A matrix  inversion 
solution  is  also  commonly  used.  Both  solutions  require  the  use  of  floating 
point  multiplication  and  addition,  while  the  forward  problem  can  be  im- 
plemented utilizing  only  fixed-point  additions.  In  this  paper,  we  develop 
a class  of  spline  basis  functions  which  solve  the  interpolation  problem  us- 
ing only  simple  arithmetic  shifts  and  fixed  point  additions  for  solutions  to 
both  the  forward  and  inverse  problems.  The  system  impulse  response  for 
the  new  interpolators  appears  to  be  closer  to  the  ideal  interpolator  than 
the  B-spline  interpolator.  We  refer  to  the  new  splines  as  (2-5-2)  splines 


§1.  Introduction 

Splines  are  well-accepted  in  Computer  Graphics[l, 2, 3, 4].  The  high  compu- 
tational requirements  of  cubic  splines  and  the  large  amounts  of  data  make 
them  difficult  to  use  in  multi-dimensional  applications.  In  many  applications, 
bilinear  interpolation  is  used  instead  of  cubic  splines  because  of  its  simplicity 
in  implementation.  However,  bilinear  interpolation  cannot  produce  images 
of  sufficient  quality  for  many  application  because  of  its  C°  continuity.  The 
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(2-5-2)  spline  approach  can  be  computed  at  nearly  the  same  efficiency  as  the 
bilinear  interpolator,  but  provides  much  higher  quality  results. 

In  general,  spline  computation  can  be  divided  into  two  parts.  The  first 
part  requires  the  solution  of  an  inverse  problem  to  obtain  the  coefficients 
(vertices)  of  the  basis  functions,  and  the  second  part  is  a forward  computation 
to  generate  the  interpolating  spline  from  its  coefficients. 

§2.  Inverse  Problem 

Assume  we  are  given  (m  — r + 3)  data  values  Pj,  % — r — 1, . . . , m + 1.  Further 
assume  that  we  wish  to  find  a continuous  spline  curve  such  that  on  a specified 

set  of  equally  spaced  points,  ux  i = r — 1 m + 1,  the  curve  attains  the 

values  Pi.  Let  Bjir(u)  denote  an  rth  order  set  of  basis  splines  defined  on  the 
knot  set  i = 0, . . . , m.  Then,  for  the  case  of  a curve  in  one-dimension  we  can 
write 

m 

Q{up)  = ViB'Atip)  = Pp,  P = k - 1, . . . , m + 1 . (1) 

t=0 

This  represents  m — r + 3 equations  inm  + 1 unknowns.  In  the  case  of  a 
cubic  spline,  r = 4 and  we  are  short  two  equations.  We  are  free  to  augment 
the  data  set  Pi  at  both  ends  with  additional  data  values.  However,  these 
augmented  data  values  effect  the  shape  of  the  interpolating  function  Q(u) 
non-locally. 

In  matrix  form,  we  obtain 

^0,r(f^r  — 2)  • • • Bm.ri'U'r  — 2) 

^0,r(tlr  — l)  Bm,r{'U'r  — l) 


B 0,r{y,m-\-2 ) •••  ^m,r(^m4-  2)-!  LKn  J L-Pm+2-J 

or 

BV_  = P, 

where  V denotes  the  vector  of  unknown  coefficients,  B represents  the  matrix  of 
basis  spline  values  at  the  knot  locations  and  P represents  the  given  set  of  data 
points  including  the  augmented  values.  The  solution  to  (2)  exists  whenever 
the  matrix  B has  an  inverse.  The  solution  to  (2)  is  efficient  whenever  B~l P 
or  its  equivalent  is  easy  to  compute. 

For  the  cubic  B-spline  defined  on  uniform  knots,  it  is  easy  to  show  that 
the  interpolation  problem  utilizing  simple,  uniform  knot  B-splines  leads  to  the 
inversion  problem 

■410.  .0 

14  10.0 

1 0 14  1.0 

6 

0.0141 
.0  . .014 


I r Pr-2  I 
Vl  Pr- 1 

• = ' (2) 
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In  this  case  the  matrix  B can  be  inverted  using  forward  elimination  and 
backwards  substitution  as  in  [1].  While  the  B-spline  matrix  is  banded  and  can 
therefore  be  inverted  reasonably  efficiently,  the  process  requires  full  floating 
point  multiplication  and  addition.  Even  after  the  matrix  inversion,  to  obtain 
the  vertices,  full  floating  point  multiplication  and  addition  must  be  used. 

The  first  attempt  at  reducing  the  number  of  the  calculations  was  made 
by  Unser,  et  al  with  a digital  filtering  solution  to  the  problem  by  noting  that 
the  inverse  of  the  FIR  filter  given  by  B(z)  = z + 4 + z~1  is  the  filter  with 
impulse  response  [3] 


b 1(n) 


—6a 
(1  - a2) 


,N 


(4) 


where,  a = \/3  — 2 is  the  smallest  root  in  absolute  value  of  the  polynomial 
z2  + 4z  + 1. 

As  they  implemented  it,  the  above  HR  filter  should  be  split  into  causal 
and  non-causal  parts,  and  applied  twice  to  obtain  the  vertices:  the  non-causal 
sequences  and  anti-causal  sequences.  Since  6_1(n)  becomes  smaller  as  |n|  gets 
larger,  we  can  assume  the  filter  has  only  several  non-zero  elements  from  the 
center  of  the  filter  sequences.  Thus,  we  can  approximate  the  filter  with  an 
FIR  filter  as  below,  and  call  it  the  Inversion  FIR  filter: 


b x(n) 


{I^a|n|,  if  M < m, 
0,  otherwise. 


(5) 


Even  with  the  Inversion  FIR  filter,  we  still  need  full  floating  point  multi- 
plication and  addition  to  get  the  vertices.  The  (1-4-1)  spline  does  not  provide 
a simple  solution  to  the  coefficient  inverse  problem  because  its  characteristic 
polynomial,  z2  + \z  + 1,  has  irrational  roots.  The  polynomials  2 z2  + 5z  + 2 has 
roots  which  are  negative  powers  of  two.  We  refer  to  these  splines  as  (2-5-2) 
splines  which  have  roots  which  are  of  similar  magnitude  to  the  roots  of  the 
(1-4-1)  spline. 

We  assume  the  spline  is  defined  on  the  uniform  set  {—2,  —1, 0, 1, 2},  and 
that  it  takes  on  the  set  of  values  {0,  |,  |,  |,  0}  at  the  knots.  In  Sect.  6,  we  set 
up  sixteen  equations  to  solve  for  the  16  polynomial  coefficients  defining  the 
(2-5-2)  spline.  We  note  that  it  has  C 2 continuity  at  knots  1,  2 and  3 and  C1 
continuity  at  the  knots  at  0 and  4 (from  the  polynomials).  Hence,  the  (2-5-2) 
spline  is  a multiple  knot  spline  defined  on  seven  knots  with  the  knots  at  0 
and  4 having  multiplicity  two.  We  also  note  that  the  normalization  property 
holds,  that  is  Xa=o  ^»(«)  = 1,  for  u € [0, 1]. 

If,  as  in  [3],  we  assume  periodic  boundary  conditions,  equation  (2)  takes 
the  following  circulant  form: 
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2 A 


A : # of  bits  in  the  fixed  point  number 
Vn:  Vertices 
Pn:  Sampled  data 

Fig.  1.  The  inverse  filter  for  a (2-5-2)  spline:  a = —1/2. 
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0 
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■V0- 
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■ Vm. 
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(6) 


or 


gB$—  ~ — ' 


(7) 


With  the  approach  IJnser,  et  al  applied  to  (1-4-1)  splines,  and  the  FIE 
filter  approximation,  we  can  obtain  a Inversion  FIR  filter  for  (2-5-2)  splines 
and  the  FIR  filter  will  be 


b-\n) 


i\M 


if  |n|  < m, 
otherwise. 


(8) 


The  filter  can  be  implemented  as  shown  in  Fig.  1 with  delays,  mulipliers 
and  adders.  Although  the  filter  shown  in  Fig.  1 can  be  used  for  a (1-4-1) 
spline  with  change  of  a = \/3  — 2,  the  FIR  filter  for  a (2-5-2)  spline  can  be 
implemented  with  an  integer  processer  because  all  its  coefficients  are  powers 
of  1/2,  and  the  power  of  two  multiplication  can  be  realized  with  shifts. 


§3.  Forward  Problem 

Ferrari,  et  al.,  provided  an  efficient  algorithm  (the  Fast  Spline  Transform, 
(FST))  for  computing  a spline  by  Derivative/Summation  [5].  Once  we  obtain 
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Z-’  Z'  Z'  Z' 

I I I I I I I I 

1st  sum  2nd  sum  3rd  sum  4th  sum 


h3 : the  3rd  order  impulse  scaling  (1,1) 
h4 : the  4th  order  impulse  scaling  (1 ,4, 1 ) 
Sn : the  interpolated  sequence 
Vn : the  vertices  sequence 
Vj’:thevertices  with  zero  padding 


Fig.  2.  The  FST  filter  for  (2-5-2)  spline  calculation. 


the  vertices,  the  spline  is  calculated  by  the  FST  with  appropriate  Impulse 
Scaling.  The  FST  is  as  follows,  where  r is  the  order  of  the  basis  function  and 
m is  the  interpolation  ratio: 

Algorithm  FST 

1)  Find  the  Dirac  functions  corresponding  to  all  orders  of  the  rth  order 
spline’s  derivatives  (r, r — 1,  • • • , 1)  [4]. 

2)  Create  an  array  of  zeros  with  m — 1 elements  between  the  knot  locations; 
initialize  k = r. 

3)  Scale  the  appropriate  Dirac  functions  by  amplitude  V„,  and  sum  the 
resulting  sequence  into  the  array  at  the  knot  locations  corresponding  to 
the  Dirac  functions(requires  shift  and  add  for  the  (2-5-2)  spline). 

4)  “Integrate”  the  resulting  sequence  once  using  repeated  summation.  Set 
k = k — 1. 

5)  Return  to  step  3.  Until  k = 0. 

6)  The  array  contains  the  values  of  the  spline  at  the  specified  locations. 

The  FST  algorithm  can  be  implemented  as  a digital  filter  for  any  spline. 
However,  to  implement  this  with  only  fixed  point  shifts  and  additions,  every 
Impulse  Scaling  element(coefficient)  of  the  spline  must  be  powers  of  two.  The 
suggested  filter  (the  FST  filter)  for  the  (2-5-2)  spline  is  shown  in  Fig.  2.  Be- 
cause the  (2-5-2)  spline  has  double  knots,  both  4th  and  3rd  order  impulse  sets 
exist,  hi  and  ha  in  Fig.  2 correspond  to  4th  order  impulse  scaling  and  3rd 
order  impulse  scaling.  Fig.  2 shows  clearly  that  the  forward  computation  of 
the  spline  is  computed  at  the  input  data  rate. 


§4.  Cosine  Examples 

Since  the  Impulse  Scaling  approach  for  the  forward  problem  always  generates 
any  (2-5-2)  spline  curves  on  the  defined  grids,  the  Impulse  Scaling  will  not 
effect  the  accuracy  of  the  generated  curves  by  each  Inverse  Problem  scheme 


Fig.  3.  The  cosine  interpolated  examples:  Top  Left:  True  Cosine  Curve,  Top 
Right:  by  Matrix  Inversion,  Bottom  Left:  by  Unser’s  Inversion  and  by 
FIR  inversion. 


(i.e.,  Matrix  Inversion,  Unser’s  HR  Filter  Inversion  and  FIR  Filter  Inversion). 
The  inversion  schemes  affect  it.  Unser’s  HR  filter  Inversion  and  FIR  filter 
Inversion  can  be  considered  as  an  approximation  operation  of  matrix  inversion. 
As  shown  in  Fig.  3,  the  cosine  curve  generated  by  matrix  inversion  is  closest  to 
the  actual  cosine  curves  because  (2)  guarantees  that  the  (2-5-2)  spline  curve 
passes  through  every  data  point  (Pn’s).  Although  the  two  filter  approaches 
are  not  guaranteed  to  pass  through  every  data  point,  they  all  produce  fairly 
accurate  curves  in  the  middle  section.  For  the  (2-5-2)  spline,  they  can  be 
implemented  with  an  integer  process  because  the  Filter  coefficients  of  HR 
filter  and  FIR  filter’s  are  powers  of  two.  It  appears  that  FIR  filter  inversion 
generates  more  accurate  cosine  than  Unser’s  IIR  filter  Invesion  in  Fig.  3,  while 
the  FIR  filter  is  nothing  but  the  approximation  of  the  IIR  filter.  The  initial 
value  estimation  for  the  anti-causal  realization  of  the  IIR  filter  results  in  less 
accuracy.  It  can  be  confirmed  that  significant  distortion  appears  at  the  right 
end  portion  of  the  IIR  filter  inversion  cosine  interpolated  curve  in  Fig.  3. 


§5.  Discussion 

With  either  FIR  Filter  Inversion  or  IIR  Filter  Inversion,  the  (2-5-2)  spline 
interpolation  will  generates  more  accurate  curves  or  surfaces  than  the  (1-4-1) 
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Fig.  4.  The  impulse  responses  of  a (1-4-1)  spline(the  dotted  line),  a (2-5-2) 
spline(the  dashed  line)  and  an  ideal  interpolator(the  solid  line). 


spline  interpolation  does.  This  can  be  easily  confirmed  by  the  fact  that  the 
(2-5-2)  spline’s  impulse  response  is  much  closer  to  the  sine  function  than  that 
of  the  (1-4-1)  spline  (refer  to  Fig.  4).  Because  the  sine  function  is  the  base  of 
the  perfect  interpolation  by  the  sampling  therom,  the  (2-5-2)  spline  produces 
more  accurate  interpolated  results  than  the  B-spline.  In  addition,  the  (2-5-2) 
can  be  implemented  with  only  fixed  point  shifts  and  additions.  Therefore,  the 
(2-5-2)  spline  is  at  least  as  accurate  as  the  B-spline,  and  can  be  computed 
much  more  efficiently. 


§6.  Polynomial  Coefficients  for  the  (2-5-2)  Spline 

The  spline  is  defined  on  five  uniformly  spaced  knots  {uo,ui,it2,U3,Ui}.  For 
each  interval,  we  define  u € [0, 1]  as  the  polynomial  variable.  Then  for  each 
interval,  the  coefficients  (a*,  hl,Ci,dl}  represent  the  polynomial  a;+feju+ CjU2+ 
diu3.  We  denote  the  four  polynomials  by  po(u),pi(u),p2(u),  and  ps(u). 

We  impose  the  following  sixteen  constraints: 


i. 

Po(0)  -- 

= 0 

ii. 

W)  -- 

= a3 

+ 

63  T C3  + 

C/4  = 

0 

iii. 

pm  -- 

= b0 

= 

0 

iv. 

pm  -- 

= 

+ 

2c3  + M3 

= 0 

V, 

Po(i)  -- 

= a o 

+ 

bo  + co  + 

dQ  = 

2 

9 

vi. 

pm  -- 

= ai 

2 

9 

vii. 

p2(i)  = 

= a 2 

■ 

fe2  + c2  + 

d-2  — 

2 

9 

viii. 

Ps( 0)  = 

= 03 

= 

2 

9 
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ix. 

x. 

xi. 

xii. 

xiii. 

xiv. 

XV. 

xvi. 

1. 

2. 

3. 

4. 

5. 


Pi(l)  — a\  + bi  + ci  + di  — | 

p2(  0)  =02  = 1 

'F’o(l)  = bo  + 2co  + 3do  = jP((0)  = fcj 
Pi'(l)  = h + 2ci  + 3di  = Pi(  0)  = b2 
P2(l)  = &2  + 2C2  + 3d2  = RJ(0)  = 63 
PS(  1)  = 2co  + 6rfo  = ^"(0)  = 2Cl 
P{'(  1)  = 2ci  + 6di  = P^'(O)  = 2c2 
P2(  1)  = 2c2  + 6 d2  = P"(0)  = 2c3 
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